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ABSTRACT 

We confirm that a high rate of shear at the center is able to stabilize an entire 
stellar disk against bar-forming (m = 2) modes, irrespective of the dark halo 
density. Our simulations of unstable power-law disks also yield the growth rate 
and pattern speed of the dominant lop-sided (m = 1) mode in close agreement 
with that predicted by linear theory in our cleanest case. The one-armed modes, 
which dominate models with extensive disks, can be controlled by tapering the 
outer disk. Two-armed instabilities of hard-centered disks are more difficult to 
identify because they are easily overwhelmed by particle noise. Nevertheless, 
we have detected the predicted m = 2 modes in simulations with very large 
numbers of particles. Such bisymmetric instabilities in these disks are provoked 
only by sharp edges and are therefore easily eliminated. 

We have constructed a cool disk model with an almost flat rotation curve 
and a quasi-exponential density profile that is unambiguously stable. The halo 
in this stable model has a large core radius, with the disk and bulge providing 
almost all the rotational support in the inner parts. 

Subject headings: galaxies: kinematics and dynamics - galaxies: halos 
- galaxies: structure - instabilities - galaxies: formation - celestial 
mechanics: stellar dynamics 
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1. Introduction 

The absence of strong bars in a significant fraction of disk galaxies has long been 
thought to present a serious problem for galactic dynamicists. Many simple models of fully 
self-gravitating disks possess global, disruptive instabilities that form bars on a dynamical 
time scale (Miller, Prendergast & Quirk 1970; Hohl 1971). In an influential paper, Ostriker 
& Peebles (1973) argued that unbarred galaxies were stable because they contained a 
large fraction of mass in a dynamically hot component which was unable to participate 
in collective instabilities. Their paper has given rise to the impression that massive dark 
matter halos are required for the stability of all disk galaxies. 

Toomre (1981) propounded a convincing mechanism for bar-forming modes. Not only 
did he provide a physical understanding of the stabilizing effect of halos for disks with 
gently rising rotation curves, but he also argued that disks could be stabilized by making 
their centers inhospitable to density waves. One way to achieve this is for the circular speed 
to remain high towards the center, which forces an inner Lindblad resonance (ILR) that 
cuts the feedback loop to the swing amplifier {e.g. Binney & Tremaine 1987 §6.3). Thus, 
even a fully self-gravitating disk can be stabilized against bisymmetric modes, much as 
Toomre's student Zang (1976) had found for the infinite, V = const, disk (here abbreviated 
as the Mestel disk, following the somewhat loose designation by Binney & Tremaine). 

Toomre's prediction was almost immediately challenged by Efstathiou, Lake & 
Negroponte (1982), who found no evidence for stabilization by hard centers in an extensive 
set of numerical simulations. The contradiction between their numerical results and 
Toomre's prediction was resolved by Sellwood (1989) who showed that large-amplitude 
disturbances in highly responsive disks could saturate the ILR, trapping particles in 
a large-scale bar similar to that which would have formed without the dense center. 
Disturbances of sufficient amplitude to trigger a bar in this non-linear manner may well have 
arisen from particle noise in the rather grainy simulations by Efstathiou et al, particularly 
when the disk is highly responsive. As a result, the stabilizing effect of a dense center has 
been often discounted, giving rise to the flawed, but frequently recited, argument that 
since maximum disks are unstable, real galaxies cannot be maximum. (A galaxy model is 
described as having a "maximum disk" when the dark matter contribution to the central 
attraction in the inner regions is negligible compared with that of the luminous disk and 
bulge.) 

Sellwood (1985, 1999) and Sellwood & Moore (1999) report A-body simulations 
of counter-examples which appeared robustly bar-stable. These realistic galaxy models 
having cool, maximal disks with dense centers, are too complicated for their stability to 
be checked by linear analysis, however. The simplest models with hard centers are disks 
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in power-law potentials, and their linear stability properties have already been calculated. 
For example, the global, small-amplitude modes in the Mestel disk were already determined 
by Zang (1976) who found, rather surprisingly, that bisymmetric (m = 2) modes were 
easily stabilized even when the disk was very cool. Zang's result also holds in other disks 
with arbitrary power-law rotation curves (Evans & Read 1998a,b). Unfortunately, all these 
power-law models suffer from strong one-armed instabilities instead, provoking Zang to 
conclude "we may merely have jumped from the frying pan into the fire." Toomre (1981) 
later pointed out that m = 1 instabilities can be controlled by reducing the disk surface 
density without changing the potential, much in the same manner as bars can be avoided 
in soft-centered disks by making them sub-maximal. 

Our objectives in this paper are two-fold. To confirm the extraordinary bisymmetric 
stability of the power-law disks, and to show that the lop-sided modes can also be quelled 
by tapering the outer disk. Such tapering again reduces the surface density everywhere, 
but changing the shape of the radial surface density profile from a power law to roughly 
exponential increases the peak disk contribution to the central attraction, preserving a 
maximum disk. 

We begin by demonstrating reasonable agreement between results from iV-body 
simulations and predictions from the linear stability analyses of power-law disks. Previous 
work has shown that the eigenmodes predicted by linear theory can be successfully 
reproduced in iV-body experiments. These tests were for the Kalnajs disks (Sellwood 1983), 
the Kuz'min- Toomre disk (Sellwood & Athanassoula 1986), and the isochrone (Earn & 
Sellwood 1995). However, all three tests were of models with uniform density cores and 
therefore finite central frequencies, and we quickly discovered that simulations of power law 
disks require much more care in order to reproduce the predicted modes. While we are 
ultimately able to claim success, our difficulties arise directly from this physical difference 
and illustrate interesting aspects of disk dynamics. 

Bisymmetric instabilities in power-law disks are provoked by a sharp inner cut-out. 
Not only are we able confirm these instabilities, but we also verify that an almost full-mass 
disk with a sufficiently gentle inner edge is bisymmetrically stable. 

We demonstrate that the lop-sided instability can be removed by tapering away the 
massive, extended outer disk. We have been able to adapt a power-law disk to construct a 
completely stable and reasonably realistic galaxy model in which the disk provides most 
of the central attraction in the inner parts. This model has a quasi-exponential disk and 
an almost flat rotation curve that arises from a combination of the massive disk, a central 
bulge, together with a halo having very low central density. We show, using both linear 
stability theory and numerical simulations, that it has no global instabilities whatsoever. 
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2. The Equilibria 



Motivated by the near-flatness of galactic rotation curves, we here study a family of 
models which possess rotation speeds that vary as a power of the radius. We break their 
scale-freeness by altering their surface density profiles using inner and outer tapers without 
changing the gravitational field. 



The circular velocity is 

2 _ ( Rq \ 2 

"^circ ~ \ R J V 0' 

where v is the velocity at the reference radius R . When (3 = 0, the rotation curve is 
completely flat; otherwise, the rotation curve is rising if (3 < and falling if f3 > 0. The 
gravitational potential in the plane of the disk is 
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Note that C{(3) = 1 for the Mestel disk, since (3 = 0. 



We will need the swing-amplification parameter X, defined by Toomre (1981) for an 
m-armed disturbance as 
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which local theory predicts is the longest unstable wavelength of axisymmetric Jeans 
instabilities in a purely rotationally supported disk (Toomre 1964). Here, k is the radial 
epicyclic frequency, as usual. 



The distribution functions (hereafter DFs) depend on both isolating integrals of 
motion: the energy per unit mass E and the specific angular momentum L z . The DFs are 

^Cr<expl-( 1 + l)E/vZ}, (3 = 



fs(E,L z ) 



(7) 
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as given by Evans (1994) for the general case and by Toomre (1977, see also Binney 
& Tremaine 1987) for the Mestel disk. In both these formulae, C is a normalization 
constant, whose value can be found in Evans (1994) or Evans & Read (1998a). The velocity 
anisotropy constant 7 prescribes the radial velocity dispersion o u of the stars: 



1 + 7 + 2f3 V R 

For these full-mass, self-similar disks, the stability parameter 



^ s 3.36 V 1+7 + 2/3' 1 } 

is independent of radius. Local theory (Toomre 1964) tells us that the disk is everywhere 
stable to axisymmetric modes provided Q s > 1. 

The DFs (^) generate the mass density implied by the gravitational potential 
through Poisson's equation. A self-similar disk with no preferred scale has unusual 
stability properties (Goodman & Evans 1999), and in any case differs greatly from the 
quasi-exponential profiles of real galactic disks. We therefore taper the disks with both 
inner and outer cut-outs, very much akin to Zang's (1976) pioneering investigation. The 
cut-out mass is still present, in the sense that it contributes to the forces experienced by 
the remaining stars, but it is not free to participate in any perturbations. The actual DF 
we use is that for the self-similar case multiplied by a taper function H(L Z ): 

f c (E,L z ) = H(L z )f s (E,L z ), (10) 

where 

TV T ft 

H(L Z ) = r- ■ ,« ir , fi ° v (11) 



L% + {v R ) v \ \L Z + Lc 

The index v controls the inner cut-out, while ft determines the outer taper; the larger 
the value of either index, the more abrupt the tapering. The outer taper is centered on a 
characteristic angular momentum L c , or radius R c , where L c = i? c t> circ (i? c ). In simulations, 
we are free to choose whatever v and jl we wish. However, the normal mode analysis is 
available only for the special values 

~ 2 + /? . 2 + f3 

V = V —p » = »—^ (12) 

where v and /1 are integers. These preferred values greatly simplify the contour integral in 
the normal mode analysis (see Appendix C of Evans & Read 1998a). Note that for the 
important (3 = case, there is no difference between v and u, ft and \x. 
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The equilibrium and linear stability properties of these models are discussed in more 
detail in Evans & Read (1998a,b). For the remainder of this paper, we adopt a system of 
units in which G = v = R = 1. Unless otherwise stated, the models in our simulations 
have v — 4, \i — 6, R c — 10 and Q s = 1. 



3. Numerical Techniques 
3.1. Aliasing 



The acceptable agreement reported in §§ |4.1| & [5| between linear theory predictions and 



our results was reached only after many failures. In the light of the successful experience 
of Earn & Sellwood (1996), we began by using the smooth-field-particle (SFP) method to 
determine forces. While less efficient than particle-mesh (PM) methods, particularly for 
larger numbers of basis functions, it does not require forces to be softened at all. 

Results with first 10 - 25 members of the &aj = 7 set of Abel-Jacobi functions (Kalnajs 
1976) were appalling. The m = 2 coefficients grew in a largely incoherent manner at a 
rate roughly ten times higher than that predicted, with "pattern speeds" which seemed 
quasi-continuous over a broad swathe of values. Non-axisymmetric features grew somewhat 
less rapidly in simulations with larger numbers of functions, but still no single "mode" 
stood out. 

After much experimentation with different time steps, particle numbers, etc., which 
yielded little in the way of improved behavior, we switched to Bessel functions as our basis 
set. Unlike the previous discrete basis, a Hankel transform has a continuous distribution of 
radial wavenumbers, k, which in practice has to be represented by a discrete set of values. 
While our first experiments revealed no significant improvement, we soon noticed that the 
average rate of growth correlated with the value of 5k adopted; i.e. the smaller 5k, the lower 
the growth rate we observed. Furthermore, inspection of the Fourier transform of the best 
fitting "mode" revealed a leading signal in excess of that predicted by linear theory, which 
gradually diminished as 5k was reduced. 

This behavior was a strong hint of aliasing. The transforms of somewhat tightly 
wrapped waves, both trailing and leading, are inadequately represented by a loose spacing 
in 5k. Thus, as noise induced waves sheared towards the trailing side, the gravitational 
field acquired components of enhanced leading waves. This purely numerical feedback 
from trailing to leading was clearly responsible for the unphysically rapid growth of 
non-axisymmetric structure in our models. The required tighter spacing in k led to rapidly 
escalating cpu time requirements. Computing significantly more than 100 functions, as our 
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results indicated would be needed, rendered this method uneconomic. (It may perhaps be 
possible to reproduce the predicted behavior with a much smaller number of functions from 
some other basis set.) 

Reverting to a PM code immediately avoided such problems, and led to much more 
physically reasonable behavior. However, particle number, grid size, time step and softening 
all needed to be substantially improved over values that seemed adequate in previous work 
before we were able to obtain results in reasonable agreement with the predictions. 

3.2. Adopted method 

We used a 2-D polar grid (Sellwood 1981) to determine the forces from the particles. 
This PM method is highly efficient but does require forces to be softened, if only slightly, in 
order to prevent relaxation (Rybicki 1972). 

We are here concerned mostly with linear instabilities, for which the radial distribution 
of mass is assumed to remain unchanged. For these cases, therefore, we need not evaluate 
the axisymmetric part of the field from the particles, and use instead the analytic central 
attraction of the infinite self-similar model. In fact, as the different sectoral harmonics of 
the perturbing forces from small amplitude disturbances are also assumed to be decoupled 
in linear theory, the particles in our linearized simulations contribute to just a single 
azimuthal Fourier component of the force field. 

As these power-law models have a large range of orbital frequencies from the center 
to the edge, it is advantageous to employ a range of time steps. The scheme adopted 
(Sellwood 1985) divides the model into a number (here 5) of radial zones, with the time 
step increasing by a factor 2 between zones. Forces arising from all the particles in one zone 
are evaluated at every step and stored separately. They are combined to find the total force 
from the entire disk, using linear interpolation between the nearest pair of times whenever 
forces are needed in the inner zones at a time not coincident with an outer step. This 
scheme is exactly time-reversible, therefore. 

As even this short step becomes inadequate near the center, we further subdivide the 
time step whenever a particle moves closer to the center than a "guard radius." For such a 
particle, we evaluate the force acting, which is anyway dominated by the fixed axisymmetric 
component this close to the center, by assuming that no others move until the next step of 
the innermost zone. The guard radius is always equal to 1/6, within which the time step is 
shorter by a factor 10. 
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The principal result from each simulation is an estimate of the eigenfrequency of the 
dominant mode. We determine these frequencies using the mode fitting procedure described 
in detail by Sellwood & Athanassoula (1986); see also Earn & Sellwood (1995). As usual, 
our estimated uncertainties indicate the entire range that encompasses all credible fits to 
subsets of the data from that run, with several scaling factors. 

3.3. Selection of Particles 

The implementation of quiet start procedures for galaxy simulations has developed 
from the initial ideas outlined by Sellwood (1983), and the latest prescription is described 
in appendix B of Debattista & Sellwood (2000). Since the DF is expressed as a function 
of the two classical integrals, f(E,L z ), we aim to create a sample of particles which has 
a density in this space as close as possible to /. We proceed by slicing accessible (E, L z ) 
space into j small areas in such a way that the integral of the DF over each area encloses a 
fraction 1/j of the total active mass. We then select (not entirely at random) a point within 
each of these areas to determine the (E, L z ) values for an orbit. For each orbit selected in 
this way, we then choose / radial phases to determine both the initial radial position and 
the azimuthal and radial velocity components of each particle. We finally place n particles 
having these three phase space coordinates at equally spaced azimuthal locations. Thus, 
the total number of particles N = jln. While as deterministic as we could make it, this 
procedure requires us to chose the shapes of each of the j areas in (E, L z )-space, and the 
values of I and n. As a result, we are able to create many different samples from the DF by 
varying these parameters. 

It seemed possible that the behavior of the modes might depend sensitively on the 
representation of the inner cut-out (see § |4.2|) . If this were true and all particles have equal 
masses, discreteness effects would ultimately become important in this region no matter 
how large the number of particles employed. We therefore experimented with packing 
more particles into this region while retaining the desired density distribution by assigning 
individual masses. In most cases, we selected particles as if there were no inner cut-out, but 
we also tried to bias the DF with an extra factor (a function either of E or of L z ) in order 
to place a still greater fraction of the particles in the inner disk. While seemingly desirable, 
this strategy actually back-fired, especially for models with declining rotation curves. By 
packing many low-mass particles in the central region, we have fewer left to represent the 
massive outer part of the disk. Results for declining rotation curves were clearly better 
when all particles had equal mass, hinting that noise in the massive part of the disk, rather 
than representation of the inner edge, is the principal source of error. 
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Parameter 



Value 



Number of particles (N) 
Grid size 



2.5 x 10 6 
200 x 256 



AR/R 
Softening length (e) 
Time step (inner zone) 



A0 



2tt/256 
~ 0.0248 
i?o/60 



0.005^0/^0 



Table 1: Standard numerical parameters in our particle-mesh simulations. 



4. Bisymmetric Modes 



4.1. Comparison with Linear Stability Theory 



Zang (1976) showed that the full-mass Mestel disk is stable to m = 2 modes unless 
the central cut-out index v > 2. For a model with v — 4 and Q s — 1, he predicted the 
eigenfrequency of the dominant m = 2 mode to be mQ p + i s = 0.879 + 0.127z; where s is 
the growth rate and fl p the pattern speed. As simulations of an infinite disk are impossible, 
we add an outer taper with index fi = 6 centered at R c = 10, which changes the predicted 
frequency very slightly to 2Q p + is = 0.879 + 0.1281 In order to keep the radial extent 
finite, we impose an additional energy truncation such that no particle's orbit takes it 
outside R = 15. This extra truncation removes little mass and is therefore unlikely to affect 
the eigenfrequency further. 

Figure [l] summarizes the results from many simulations of this one model as the 
numerical parameters are varied. Each pair of panels shows how the measured pattern 
speed and growth rate vary as one numerical parameter is changed from the "standard" 
values given in Table [I]. The horizontal dotted lines indicate the frequency predicted from 
linear theory. 

As the number of particles is increased (panels a), the measured frequencies clearly 
approach the predicted value near N = 2.5 x 10 6 . The scatter results from making many 
different selections of particles from the DF for each value of N. The steady trends seen in 
panels (b), (c) & (d), on the other hand, result from models which all began with exactly 
the same sample of particles. 

The data we fit to estimate the eigenfrequency are never well approximated by a single 
mode, and all our fits require a second "mode." Unlike the dominant mode, however, 
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Fig. 1. — The effects of varying the numerical parameters on the pattern speed and growth 
rate of the dominant bisymmetric mode of the Mestel disk with v = 4, // = 6 and R c = 10. 
Our standard values are given in Table |l|. The parameters that are varied are: (a) number 
of particles N (no error bars shown), (b) softening length e, (c) grid resolution A0 and (d) 
time step St. The prediction from linear theory is shown as a dotted line. 
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Fig. 2. — The effects of varying the inner cut-out index v on the pattern speed and growth 
rate of the dominant bisymmetric mode of the Mestel disk. The circles show the linear theory 
predictions and the error bars show the dispersion from several simulations with different 
samples drawn from the DF for each case. 



the second wave does not have a coherent frequency throughout its growth, neither is its 
frequency similar from runs with different samples from the DF. Since amplified particle 
noise can be reasonably well approximated by a coherent wave for short time periods, 
our fitting procedure is never able to separate noise completely from the true mode. The 
scatter in panels (a) and the lack of scatter in the other panels strongly suggest that our 
measurements are most affected by the noise spectrum in the simulation, which is changed 
with every new sample drawn from the DF. We discuss particle noise in more detail in §IO 



The trends all show that the results have yet to "converge" as each parameter is 
refined. Our standard parameters (Table |l|) are not perfect, but seem adequate. It may 
seem worrisome that the trends in panels (b), (c) & (d) are not converging to the prediction. 
We stress, however, that these results were obtained from a single sample of particles, and 
the true error bars on these results should be enlarged considerably to take account of the 
~ 5 — 10% scatter for this N in panels (a). 

The scatter in panels (a) seems to decrease as N rises suggesting that the predicted 
values might be reliably reproduced in still larger calculations, as they should. We obviously 
should expect our higher quality results to cluster around the prediction, which we claim 
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Predicted Frequency 
2Q p + is 


Observed Observed 
Real part Imaginary part 


n ok 
-u.zo 


u.ooo t" u. ry it 


U.o^tU It U.Ulo U.l / 1 I u.uzu 

0.831 ±0.001 0.186 ±0.003 
0.832 ±0.004 0.189 ±0.010 
0.864 ±0.001 0.198 ±0.005 


0.25 


0.920 + 0.060i 


0.757 ±0.011 0.097 ±0.017 
0.865 ±0.038 0.085 ±0.011 
0.805 ±0.013 0.102 ±0.010 
0.865 ±0.002 0.086 ±0.004 



Table 2: Predicted and estimated eigenfrequencies of m = 2 modes in power-law disks. The 
upper panel refers to the models with rising rotation curves (/3 = —0.25), the lower panel to 
falling rotation curves (/3 = 0.25). 

they do. 

The predictions and results for several values of the inner cut-out index v are shown 
in Figure |2|. The error bars on the data points show the 1-a dispersion (not the standard 
error) from a number (> 4) samples from the DF in each case. All these simulations 
used the standard numerical parameters given in Table [I]. We find reasonable agreement 
between linear theory and the simulation results, but generally far below the level Earn & 
Sellwood (1995), who employed one twentieth the number of particles, were able to achieve 
for the isochrone disk. Even with the larger N used here, no bisymmetric mode, however 
strongly growing, stands out clearly from the noise in the manner observed routinely in 
simulations of soft-centered disks. The values we obtain have a tendency to underestimate 
the prediction, of Q p in particular, which we find bothersome; we have been unable to 
eliminate this small systematic difference and believe residual noise to be the most likely 
culprit. 

Our results for two different power-law models, one with rising and the other with 
falling rotation curves, are summarized in Table together with the predictions from 
linear theory, for v = 4 in each case. The agreement for the model with the rising rotation 
curve (j3 = —0.25) is satisfactory, but that for the falling rotation curve model is much 
poorer. Further investigation seems to indicate that the number of particles needed for 
good agreement rises with increasing (3. One part of the reason for this behavior was traced 
to the folly of filling the inner cut-out with a large fraction of low-mass particles, as noted in 
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§ |3.3| . But this is clearly not the full story, since the results given in the Table |^ (for this case 
only) are for equal mass particles. While slowly-growing modes are clearly more difficult to 
reproduce, § [4.3| gives a further reason why the falling rotation curve case is worse. 

4.2. Mode mechanism 

The reason the modes of soft-centered disks are so much easier to reproduce than those 
of the power-law disks stems from the difference in reflection mechanisms between the two 
types of disk. In both cases, they are swing-amplified inner cavity modes (Toomre 1981) in 
which in-going trailing waves return as out-going leading waves providing positive feedback 
to the amplifier. A wave-train in a soft-centered disk, such as the isochrone, reflects off 
the disk center and returns as a leading wave without attenuation. The centers of the 
power-law disks, on the other hand, cannot support density waves for two reasons: first, 
frequencies become arbitrarily large as R — > 0, which ensures an ILR at some radius for all 
pattern speeds, and second, the inner cut-out causes Q to rise steeply towards the center. 
Feedback in this disk is through reflection off the inner cut-out which probably returns some 
signal in both the short-leading and long-trailing wave channels. It is also possible that 
some fraction of the in-coming wave is absorbed by the imperfectly shielded ILR, especially 
when the resonance is broadened by the rapid growth of the mode. The sharper the inner 
cut-out, the stronger the reflected leading wave becomes and the more vigorous the mode. 

Thus, unlike for the soft-centered disks, only a fraction of the in-going trailing wave 
returns as out-going leading wave to provide feedback to the swing amplifier. The fraction 
returning as long-trailing waves also shears to short-trailing near corotation (e.g. Mark 
1976), but with mild amplification only when Q < 1.2. The effect of this "leakage" from 
the main swing- amplified feedback loop is to slow the rate of growth of the mode. The 
response of the disk to noise, on the other hand, is unaffected by feedback and is therefore 
proportionately more troublesome in the power-law disks. 

4.3. Amplified particle noise 

When the reflected leading signal is weak, as in disks with gentle cut-outs, the predicted 
mode is easily submerged by particle noise, but particle noise remains an unusually severe 
problem even in disks with sharp cut-outs. The leading component of the noise spectrum 
produces apparently growing non-axisymmetric features that obscure our view of the 
underlying true instability. 
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Fig. 3. — Left: Radial variation of the fraction of mass per square A cr i t , £, for the three 
different power-laws: the dashed curve is for (3 = 0.25, unbroken for (3 = and dotted for 
(3 = —0.25. The number of particles per A^ rit is N times the quantity plotted. Right: The 
radial variation of Q for these three disk models. 

When particles are distributed randomly, the level of particle noise depends both on 
the number of particles per square A cr i t {e.g. Toomre & Kalnajs 1991) and the vigor of 
the swing amplifier. The maximum growth that the swing amplifier can produce depends 
on several properties of the disk: the parameters Q and X, as well as the shear rate 
r = 1 + (3/2. We find that the expected amplification remains roughly constant as (3 varies. 
(We held Q s constant and the effects of changing both F and X largely cancel.) 

The graininess of the disk, on the other hand, does depend on the power-law index, as 
shown in Figure |^. When all particles have equal mass, the number of particles per square 
A cr it is N£ where £ is the ratio of mass per A^ rit to the total active mass in the simulation. 
The regions of the disk where £ is small have high Q (as shown in the right-hand panel) and 
do not therefore cause difficulties. Note that the peaks of these curves are of order unity, so 
we have millions of particles per A^ rit in our largest calculations. 

A quiet start is designed to reduce the seed amplitude of non-axisymmetric features 
in order to be able to observe linear growth over many e-folds. By placing a number of 
particles at equal angular spacing on rings and restricting forces to low order sectoral 
harmonics, initial disturbance forces are many orders of magnitude below those that would 
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Fig. 4. — Snapshots showing the position of one particle in 500 from a noisy start simulation 
with N = 2.5 million particles in which perturbation forces were restricted to m < 4 and 
m/1. The radius of the circles is 17.2R and the times are in units of Rq/vq. The model has 
no true global instabilities. Swing-amplified noise quickly produces large-amplitude patterns 
and later a bar even for this large N. 
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have arisen from random initial positions. Such an arrangement is not entirely noise free, 
however. We provide an initial seed noise spectrum by nudging every particle away from 
perfect azimuthal spacing by a randomly chosen fraction of 0.002 radians, or about 0.1°. 
This initial noise seeds the growth both of the desired global mode (if one is present) 
as well as a swing-amplified response. The resulting enhanced disturbance forces further 
distort the regular particle arrangement, creating enhanced density fluctuations which are 
further amplified, and so on. Thus the quiet start is itself unstable, and non-axisymmetric 
amplitudes rise until the distribution of particles is essentially random. The rate at which 
swing-amplification drives the break up again depends on the responsiveness of the disk as 
well as iV£. 

In order to quantify the rate at which quiet starts disrupt, we have run simulations 
of our usual Mestel disk but with the inner cut-out index v — 2. Zang (1976) concluded 
that the Mestel disk is well stable to modes with m ^ 1 in this case and our addition of an 
outer taper does not alter his finding. Swing-amplification is still present in this "stable" 
case, however. The amplitude of m = 2 disturbances in a quiet start simulation grew in a 
quasi-exponential manner at an approximate rate of 0.06, but we could find no evidence for 
coherent modes. This simulation confirms that swing-amplification from the initial noise 
causes a rapid break-up of the quiet start even when the equivalent smooth disk has no 
instability. 

Amplified particle noise is the reason that growth rates appear to increase as iV 
decreases in Figure |l](a). Only when N > 10 6 is the break-up rate slow enough to have 
little effect on the growth rate of the underlying mildly growing normal mode. For fixed N, 
the P = 0.25 case has the smallest number of particles per A^ rit (dashed curve in Figure 0), 
and amplified particle noise should be the greatest nuisance for this case, as we have found. 
It should be noted that quiet starts neither create nor alleviate this difficulty because they 
reduce the amplitudes of both the mode and the noise in proportion. 

4.4. Confirmation of bisymmetric stability 

The severity of the noise problem is illustrated in Figure [| which shows snapshots from 
a noisy start simulation of the linearly stable model; numerical parameters were as given 
in Table [TJ except that we placed each of the 2.5 million particles at a randomly chosen 
azimuth and included forces from the density sectoral harmonics through 4, with only 
m = 1 omitted. Despite there being no true global instability, after just 40 dynamical times 
(or one orbital period at R ~ 6), swing amplified particle noise has already produced a 
visible spiral pattern of sufficient amplitude to alter the density distribution and heat the 
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inner disk. 

Noise-driven activity in this lively Q s = 1 disk is so strong, even for this large number 
of particles, that the model later develops a strong bar. The bar does not result from a 
linear instability, but through saturation of the ILR, or non-linear trapping, as discussed by 
Sellwood (1989). The swing-amplifier is tamed by raising Q s and a further long simulation 
has confirmed the global stability of a similar model with Q s = 1.5. 

5. One-armed modes 

5.1. Comparison with Linear Stability Theory 

Linear theory (Zang 1976; Evans & Read 1998b) predicts a dominant pair of m = 1 
modes with closely spaced frequencies for cool power-law disks. For our centrally cut-out 
and outer truncated Q B — 1 disks, the predicted frequencies for the dominant mode pair are 
given in Table |3|. These frequency differences are typical. The Table also presents results 
from a number of separate simulations using different particle samples, but were otherwise 
identical for each disk and used the numerical parameters given in Table |l|. 

In all our simulations of the Mestel and power-law disks, the m = 1 data were consistent 
with just a single growing mode having the indicated eigenfrequency, but we were able to 
find hints of the two predicted modes in one case. Separation of two such closely-spaced 
modes would require a long period of integration in the simulations, which cannot be 
continued indefinitely as the modes saturate after approximately 8 e-folds, or roughly 100 
dynamical times. In the linear regime, therefore, they will have rotated relative to one 
another by only some 1.6 radians. It seems reasonable that we should find growth rates 
consistent with that of the more rapidly growing mode, but it is curious that we sometimes 
find a pattern speed closer to that of the more slowly growing mode. The variations within 
a set of experiments using the same physical model probably arise from different relative 
amplitudes of the two modes seeded by the different initial particle noise. 

5.2. Discussion 

The reason for the doublet of m = 1 modes is that there are two possible paths from 
trailing to leading waves in the inner disk (Toomre 1998, private communication). This 
is because the reflection of an in-going trailing wave off the inner cut-out returns some 
short-leading and long-trailing signal. Both paths complete the feedback loop through 
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R 
P 


Predicted frequencies 

MnHp 1 MnHp 9 


Observed Observed 

rvcdi udi i midciiiidi y udi i 


-0.25 


0.176 ±0.130i 0.140 + 0.089i 


0.163 ±0.004 0.128 ±0.004 
0.162 ±0.003 0.125 ±0.002 





0.189 + 0.118i 0.173 ±0.089i 


0.172 ±0.002 0.105 ±0.002 
0.169 ±0.003 0.111 ±0.007 
0.173 ±0.005 0.118 ±0.006 


0.25 


0.213 + 0.100i 0.198 + 0.081i 


0.200 ±0.006 0.098 ±0.005 
0.200 ±0.002 0.095 ±0.001 



Table 3: Predicted and estimated eigenfrequencies of m — 1 modes for disks with three 
different power-law indices (5. 




Fig. 5. — The points with error bars show estimates from simulations of (a) the pattern 
speed and (b) the growth rate of the dominant lop-sided mode in a Q s = 1.5 Mestel disk. 
The prediction from linear theory is shown by the dotted lines. 



swing amplification at corotation. As Q increases, the distinction between the two paths 
diminishes, and the splitting of the doublet decreases. When Q ~ 1.5, linear theory recovers 
a single m = 1 mode. 



In order to verify that the discrepancies between predictions and simulations are 
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due to the existence of the closely spaced doublet, we have tested the prediction of an 
isolated m = 1 mode in a Q s = 1.5 Mestel disk. Unlike all results reported above, in this 
case the data from a simulation are well approximated by a single coherent mode, with 
little interference from noise. Furthermore, we no longer find significant variations in the 
measured frequency between different draws of particles when N = 2.5 x 10 6 . Figure ||] 
shows that in this clean case we are able to achieve truly convincing agreement between 
linear theory and simulation - the discrepancy is less than 2% for the largest iV we have 
used. This result suggests that in this case almost all the in-going trailing wave reflects to 
a leading wave, thereby allowing the mode to outgrow any interference from particle noise 
and to stand out quickly and cleanly. The cleaner reflection at m = 1 is probably largely 
due to the greater spatial scale of the mode, but the absence of an ILR for m = 1 waves 
may also be significant, since no residual damping is possible. 

6. A Completely Stable Disk Model 

6.1. Linear Stability Theory 

Our aim here is to provide a clear-cut example of a stable, massive disk with a nearly 
exponential radial profile. We adopt a power-law disk with a gently falling rotation curve 
(P = 0.05), but choose the inner cut-out index u = 1 and the outer taper index \i = 2. 
Figure || shows how the pattern speed and the growth rate of the dominant one-armed mode 
depends on the characteristic radius of the outer taper R c for a disk of nominal Q s = 1.0. 
[Note that the nominal Q s value given in equation (||) refers to the underlying scale-free 
disk, whereas the actual Q value of the tapered disk is higher.] Figure ^| shows that when 
the outer taper is centered on R c < 5, this disk becomes completely stable to one-armed 
modes. We have verified that this disk is also linearly stable to the m = 2, 3 and 4 modes. 
At least within the framework of linear theory, the strongly tapered disk is completely 
stable, therefore. 

6.2. A quasi-exponential disk 

We have removed material from the disk with inner and outer tapers without changing 
the potential. The extra central attraction is provided by unresponsive matter, which we 
think of as two rigid spherical components, a bulge and halo, whose density distributions 
are determined as follows. 
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Fig. 6. — The pattern speed Q p (dashed line) and growth rate s (full line) of the dominant 
m = 1 mode as the position of the outer taper R c is varied. In this Q s = 1.0 power-law disk 
with f3 = 0.05 ,v — 1 and /i — 2, the model is stable against one-armed instabilities when 
Rc, < 5. 



The disk surface density has a logarithmic spiral decomposition: 
T,(R) = S / daB(a) exp (ia - |) log 

J— CO ' 

where B(a) is given by 
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For the power-law disks, the integration is analytic; the details are given in Appendix A. 
With this decomposition in hand, we find the spherical density pd that reproduces the 
forces in the equatorial plane provided by the disk is 

„ 5 



z_in / T \ ~ 2 f°° 

P ^ r) = W \Rj Jo da (a 2 + l)K(a,0)B(a)exp 



ia log ( — 

Kq 



(15) 



where K(a, 0) is the axisymmetric Kalnajs (1971) gravity function. Therefore, the spherical 
bulge density pb and the halo density ph must satisfy 



Pb(r) +Ph(r) 



-Pd[r), 



(16) 
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Fig. 7. — (a) The rotation curves of the disk (unbroken line), bulge (long-dashed) and halo 
(short-dashed) are plotted against cylindrical radius R. The total rotation curve (dot-dashed) 
is gently falling. Throughout the inner galaxy, the disk and bulge are the major contributors 
to the rotation curve, (b) The radial variation of the disk surface density (unbroken line) 
which is quasi-exponential. The dot-dashed line shows the surface density of the Mestel disk 
for comparison, (c) The true Q profile of the disk (unbroken line) is shown against radius 
R. It rises sharply near the center and in the outer parts, although there is an intermediate 
regime where Q ~ 1.75. The dot-dashed line shows the nominal Q s — 1.0. (d) The density 
of the spherical bulge (long-dashed line) and halo (short-dashed line) are shown against 
spherical radius r. 
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with pd on the right hand side given by equation (15). There is some freedom in the 
partition of the rigid density between the bulge and the halo. We choose the dark halo to 
be a cored power-law sphere with the potential-density pair (e.g. Evans 1993) 

«V(l-/3)+3<i) = vj m 

To ensure the halo makes a modest contribution to the rotation curve in the inner galaxy, 
we chose the core radius to be large, = 20, giving it an almost uniform density in the 
inner parts. At large radii, on the other hand, the halo provides almost all of the rotational 
support. The remainder of the cut-out density is ascribed to a bulge. 

The equilibrium properties of this model are shown in Figure |7|. Panel (a) shows the 
rotation curves of the disk, bulge and halo components. The whole model has a slowly 
falling rotation curve. The active disk has a quasi-exponential profile (panel b). The 
distribution functions given by equation (|10|) may be the simplest known which yield 
quasi-exponential disks with flat rotation curves. As the true Q profile of the disk (panel 
c) is always larger than unity, the disk is axisymmetrically stable. At both small and large 
radii, Q rises sharply where a small fraction of the mass remains active, while Q < 2 at 
intermediate radii. Panel (d) shows the density of the spherical bulge and halo. As the halo 
makes a negligible contribution to the mass budget in the inner parts, we can regard this is 
a "maximum disk" model, albeit with a massive bulge. 

It should be noted that although the surface density (Figure 0b) reaches at most barely 
more than half that of the untapered disk near R ~ 2, the tapered disk's contribution to the 
circular speed (Figure [F|a) peaks at some 75% of that from the full-mass disk near R ~ 5. 
This result should not be surprising since the taper has cut away material at larger radii 
which previously depressed the disk's contribution to the central attraction. 



6.3. Simulations 

We have used our iV-body code to test the linear theory prediction that this model 
has no global instabilities. Figure || shows the result from a noisy start simulation with a 
modest number of particles (N = 5 x 10 4 ) and with m = — 4 forces calculated from the 
particles. All other numerical parameters were as given in Table The force acting on each 
particle at every step included the fixed central attraction of our bulge/halo in addition to 
the forces computed from the particles. 

The simulation quickly develops some spiral patterns, which result from swing-amplified 
particle noise, but there is no evidence for a coherent inner bar or any lop-sidedness. A test 
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Fig. 8. — The positions of one particle in 10 from a noisy start simulation of model described 
in § |6.2| . The radius of the circles is 17.2-Rq and the times are in units of Rq/vq. The simulation 
employed N = 5 x 10 4 and forces from m = — 4 were all active. The model appears to 
be stable to bar-forming and lop-sided modes, but as usual swing-amplified noise quickly 
produces spiral patterns which heat the disk, especially in the inner parts. 
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with ten times more particles behaved similarly, but the spirals were so weak that the disk 
barely heated at all, and no long-lived coherent waves were detectable at any m. 

We conclude that this model is robustly stable. 

7. Conclusions 

It should no longer be doubted that real galaxies can be stabilized by dense centers, 
as argued by Toomre (1981). We have finally confirmed Zang's (1976) global stability 
analysis for the Mestel disk, and also its generalization to other power laws by Evans & 
Read (1998b). This encouraging agreement between linear theory and simulations for these 
idealized galaxy models buttresses Sellwood's (1985, 1999) claims of global stability for 
more realistic, but less analytically tractable, models with hard centers. 

The lop-sided (m — 1) instabilities of the power- laws disks, first discovered by Zang, 
can be controlled by reducing the disk surface density to < 2/3 of that of a full-mass disk, 
without changing the potential (Toomre 1981), much in the same manner as bars can be 
avoided in soft-centered disks by making them sub-maximal. In this paper, we have been 
able to show that the m = 1 mode can also be avoided by tapering the massive outer 
disk until its surface density declines in a quasi-exponential manner. The mechanism is no 
different since the surface density reaches barely more than half that of the untapered disk 
at any radius, but the result is a more realistic model in which the disk contribution to the 
central attraction in the inner parts is decreased to a smaller extent than is the mass. 

Many real galaxies do appear to be stabilized by dense centers. High-resolution 
kinematic observations of spiral galaxies by both Rubin, Kenney & Young (1997) and by 
Sofue et al. (1999) have revealed high orbital speeds close to their centers. In the cases 
reported by Sofue et al, for example, the only galaxies with gently rising rotation curves 
all have rather low luminosity (a peak circular velocity < 150 km s" 1 ). Such galaxies are 
required to have large dark matter fractions (e.g. Broeils 1992) - presumably enough to 
suppress bar-forming instabilities. Almost every galaxy with a circular speed in excess of 
150 km s -1 has a steep inner rise in the rotation curve, and must therefore be bar-stable 
whatever its dark matter content. 

The linear modes of these power-law disks proved to be much harder to identify than 
those in previous studies. The predicted bisymmetric instabilities are completely obscured 
by swing-amplified noise unless we employ much larger numbers of particles than are needed 
in soft-centered disks, and even then the underlying mode does not stand out clearly. We 
attribute this difference to the different reflection mechanism in the feed-back loop for 
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the power-law disks. The mode stands out clearly from the noise only in cases where the 
inner reflection is total, as in soft-centered disks, and possibly also for m — 1 modes in 
hard-centered disks (§5.2). But particle noise readily swamps the weakened modes which 
result when the reflection involves considerable leakage, as happens for the bisymmetric 
modes in hard-centered disks (§4.2). 

We have been able to confirm the predicted most unstable m — 1 mode in the Q s = 1.5 
disk with a precision approaching that obtained by Earn & Sellwood (1995) for the bar 
modes in the isochrone disk. Our simulations of the cooler Q s = 1 disk appeared consistent 
with only a single mode where linear theory predicted a closely-spaced doublet; nevertheless, 
the frequency we obtained was always some "average" of the predicted pair. 

We have shown that one-armed instabilities are tamed by tapering the outer disk. We 
note that the lop-sided instability has implications for the dark matter candidate proposed 
by Pfenniger, Combes & Martinet (1994). They suggested that flat rotation curves result 
from a high surface density of very cold molecular gas in a rotationally supported disk 
extending beyond the optical edge. This hypothesis is ruled out on stability grounds: 
massive extended disks with flat rotation curves are grossly unstable to m = 1 modes. 

We have presented a realistic galaxy model that is globally stable. It has a quasi- 
exponential disk in an almost flat rotation curve potential, and includes both a halo with 
a large core radius and a bulge. The model is stabilized against bar modes by its dense 
center and against one-armed modes by the exponentially decreasing surface density profile. 
Even though the halo makes an negligible contribution to the central attraction in the 
inner parts, the model is clearly globally stable as we demonstrate using both normal mode 
analysis and iV-body simulations. 

We are indebted to Alar Toomre for his comments on an early draft of this paper and 
detailed, patient explanations of the mode mechanisms. Useful conversations with David 
Earn are also acknowledged. We thank the referee (Agris Kalnajs) for a thoughtful and 
stimulating report. This work was begun at the Isaac Newton Institute, Cambridge and 
continued at the Lorentz Centre, Leiden. The hospitality of both places is acknowledged. 
JAS is supported by NSF grant AST 96/17088 and NASA LTSA grant NAG 5-6037. NWE 
thanks the Royal Society for financial support. 
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In section |6.2| , we exploit the decomposition into logarithmic spirals to find the 
three-dimensional spherical bulge density that provides the forces to compensate for 
the matter carved out of our disks. The integration for the logarithmic spiral density 
transform ([14]) can be carried through explicitly for the power-law disks. This is done by 
writing the density as an integral over velocity space of the cut-out distribution function, 
and then reversing the orders of integrations. 

For the case of our quasi-exponential disk, the inner cut-out is v — 1 while the outer 
taper is \x = 2. The transform then is 



A. The Density Transform 
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This is valid for > 0; slightly different expressions hold good for < 0. 
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